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Abstract: We derive a new method for initial-state collinear showering in Monte-Carlo event 
generators which is based on the use of unintegrated parton correlation functions. Combined 
with a previously derived method for final-state showering, the method solves the problem 
of treating both the hard scattering and the evolution kernels to be used in arbitrarily non- 
leading order. Although we only treat collinear showering, so that further extensions are 
needed for QCD, we have discovered several new results: (1) It is better to generate exact 
parton kinematics in the hard scattering rather than with the subsequent parton showering, 
and similarly at each step of the showering. (2) Parton showering is then done conditionally 
on the exact energy-momentum of the initiating parton. (3) We obtain a factorization for 
structure functions in terms of parton correlation functions so that parton kinematics can be 
treated exactly from the beginning. (4) We obtain two factorization properties for parton 
correlation functions, one in terms of ordinary parton densities and one, suitable for event 
generation, in terms of parton correlation functions themselves. 
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1. Introduction 



Monte-Carlo event generators (MCEGs) provide a very powerful approach to implement- 
ing perturbative QCD to predict many observables in hard processes. However, despite a 
considerable amount of work — see the paper by Frixione and Webber Q and the refer- 
ences reviewed there — current algorithms do not systematically extend beyond an improved 
leading-logarithm approximation. This is in contrast to factorization theorems for hard inclu- 
sive processes, where it is known how both hard-scattering coefficients and evolution kernels 
may be computed, in principle, to arbitrarily high order in perturbation theory. As a compu- 
tational tool, event generators not only provide estimates of the exclusive components of cross 
sections, but they have the important property that event generation needs a time roughly 
linear in the number of particles. In contrast, the computation of Feynman graphs for cross 
sections involving N particles needs resources proportional to roughly N\. 

In this paper, we provide a solution to the problem of non-leading corrections in MCEGs 
when initial-state collinear showering is involved. Combined with a method one of us devel- 
oped |2j for final-state collinear showering, this allows a complete treatment for deep-inelastic 
scattering, for example. Our treatment in this paper is restricted to collinear showering, so it 
applies as it stands only to non-gauge theories. However, even though substantial extension 
of the method will be needed to handle QCD, with its soft gluon emission, our work solves a 
number of important conceptual problems whose solution is also needed in QCD. 

We will give a more detailed account of the problems we are concerned about in Sec. 
^. Here we will simply remark that the algorithms in standard MCEGs are based on the 
factorization theorem for inclusive hard processes, whereas MCEGs concern themselves 
with the exclusive components of cross sections. Thus the standard factorization theorem does 
not really provide a fully satisfactory foundation for MCEGs; a more general and powerful 
theorem is needed. This can be seen in several ways: 

• The proofs ^ of standard factorization explicitly involve sums over unobserved final 
states, but MCEGs resolve the full final state. 

• Approximations on the kinematics of internal partons are made that become invalid 
if the exact final state is observed: the approximations violate exact momentum con- 
servation. Moreover, as we will show in Sec. |2|, the approximations interfere with the 
factorized structure of probabilities (or weights) that enables the stochastic generation 
of events to take a time linear in event size. 

• In standard factorization, cross sections that are infra-red and collinear safe are approx- 
imated by singular distributions like 

(1 + a s A) 5(x) + a s B (^j +.... (1.1) 

While such approximations can be valid when integrated with a smooth function, they 
are not valid when a fully differential cross section is under consideration. 
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• Although the previous problem can be alleviated by some kind of resummation, it is 
preferable to have an improved factorization theorem from the beginning. 

The most fundamental part of the proof of factorization, the analysis of the momentum-space 
regions that give leading contributions to the cross section, is indeed generally applicable, 
and so it does give a correct starting point for deriving algorithms for MCEGs. What is at 
issue are the further steps, and in particular the approximations, that are used to obtain the 
standard form of factorization; these approximations are only appropriate when there is an 
inclusive sum over the hadronic final state. 

Our approach solves the problems by working with parton correlation functions that are 
defined without an integral over parton kinematics, unlike conventional parton densities. This 
avoids the need to use approximated parton kinematics in places where the approximation 
is invalid. Moreover, we will also work with parton correlation functions and cross sections 
that are defined without an integral over hadronic final states (or, more conveniently, with 
an arbitrary weight function of the hadronic final state). This gives us the appropriate tools 
for systematically and precisely treating the actual physical situations which MCEGs aim to 
approximate. 

Of course, parton correlation functions are more complicated objects than parton densi- 
ties, so without further information a lot of predictive power is lost. We solve this issue with 
extra factorization theorems for the parton correlation functions at large parton virtuality; 
the showering algorithm is what implements these. 

At low parton virtuality, the information in parton correlation functions corresponds to 
the information in non-perturbative hadronization models that are present in conventional 
MCEGs. There is no gain or loss of predictive power here. 

We will find that we can explicitly represent all probabilities and cross sections for a 
MCEG in terms of integrals over parton correlation functions, and hence in terms of explicit 
operator matrix elements. The systematic extension to non-leading order and non-leading 
logarithms then becomes straightforward, with the use of the subtractive methods advocated 
by us in g |. 

Related work is in papers by Watt, Martin and Ryskin || ||| who have proposed the use of 
what they call "doubly unintegrated parton distributions" (DUPDFs). This concept coincides 
with our "parton correlation functions". Their primary motivation is the same as ours: the 
need to treat parton kinematics exactly in certain observables in a way that is fully compatible 
with ordinary parton densities. Their work treats QCD, so it is of direct phenomenological 
significance. However, it is also restricted to a kind of leading approximation; they obtain a 
result for the DUPDFs from an examination of the last step of the evolution. This corresponds 
to one of our extra factorization theorems for parton correlation functions, in its leading- 
order approximation when the parton kinematics are strongly ordered. Our aims are more 
ambitious even if we have not achieved them in QCD so far: we have obtained a formalism 
that can be applied to arbitrarily non-leading corrections and in the context of a MCEG. 
To this end, we have provided exact operator definitions of our distributions. In QCD, as 
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we will explain, it is necessary to deal with the issues symptomized by the gauge invariance 
problem of the definition of the pCfs: Where should the Wilson lines be? In a leading 
logarithm approximation, these issues are readily evaded by inserting appropriate cutoffs, 
without examining exactly how these are to be implemented in an exact operator definition. 

The organization of the paper is the following. After explaining the conceptual issues 
that concern us in Sec. 0, we describe our model theory in Sec. [| Then in Sec. || we derive 
our new factorization for DIS structure functions in terms of parton correlation functions. In 
Sec. ||, we obtain factorization theorems for parton correlation functions; these give both the 
showering algorithm and the computation of parton correlation functions in terms of ordinary 
parton densities (as in the MS scheme). In Sec. ^ we show an algorithm for a MCEG that 
implements our factorization formulas. Then in Sec. [7| we give explicit calculations of the 
next-to-leading order (NLO) hard scattering for deep-inelastic scattering. Finally, in Sec. ||, 
we give our conclusions and an indication of the approach we propose to use to extend our 
results to full QCD. 



p 

Figure 1: Generation and structure of event from MCEG for deep-inelastic scattering. One vertex 
has four attached lines; it would arise from a possible NLO branching. 

2. Rationale for a new algorithm 

The overall structure of a MCEG is that each event is generated by the following steps: 

1. Generate a hard scattering. This results in a set of partons and values of their momenta 
(or of some scalar variables that can be used to determine the momenta). 

2. For each available parton, either choose not to shower it, or generate a branching, i.e., 
a set of new partons. 

3. Go back to the previous step, applying it to any new partons that have been generated. 

This general structure corresponds to a factorization theorem, and the generated events can 
be given a Feynman-graph-like structure, as in Fig. [I] for deep-inelastic scattering (DIS). 
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The generation of the hard scattering is represented by the vertex labeled q, and the resulting 
partons are represented by the attached lines. The branchings of these partons are represented 
by the vertices at the other ends of the lines. The arrows label momentum flow with respect 
to a corresponding Feynman graph vertex. The lines with no vertex at the outgoing end 
are final-state particles, and the one line, labeled P, with no vertex at the incoming end is 
the target particle. Some generalization is needed in the case of angular ordered coherent 
branching, but this does not affect the overall ideas. Notice that the lines represent both the 
generated partons and the flow of information between different units of computation, at the 
vertices. 

The key property that gives a straightforward 
algorithm is that at each step the showering proba- 
bility for generating partons depends on a fixed num- 
ber of variables that correspond to the single parent 
parton of that step. The showering probability is, of 
course, a conditional probability. 1 Since at each step 
we have several partons, each of which needs to be 
showered, the showering probabilities for the differ- 
ent partons should be independent; this corresponds 
to a factorization theorem. What is needed to allow 

a clear systematic discussion of non-leading correc- Figure 2: Deep inelastic scattering, 
tions is a definite operator matrix element expres- 
sion for the probability for showering; any non-trivial communication between the showers, 
as exhibited by a non-factorization, will tend to preclude the possibility of a matrix element 
formulation. Without a systematic factorization, 2 so that different showers are independent, 
it is hard to give a systematic treatment for the production of arbitrarily many particles. 

As we will now show, the most obvious formulation for the algorithm, founded on the 
parton model approximation to factorization, actually violates some of these properties be- 
cause of its treatment of parton kinematics. This is what motivates us to find a more general 
approach. 

For ease of presentation, we shall work with a model theory c/> 3 in space-time dimension 
n = 6 throughout the paper, the methods immediately generalize to any nongauge theory. 
The details of the model will be described in Sec. |3[ 

We consider DIS, whose kinematics are shown in Fig. ^. There q is the photon momentum 
and P the hadron momentum. We work in the virtual-photon-hadron center-of-mass frame, 
using light-cone coordinates, = (p ° ±p z )/y/2, where, up to corrections power-suppressed 

1 In a subtractive formalism beyond leading order, such as we will use, correction terms to the cross section 
may sometimes be negative, so we will allow the generation of weighted events, some with negative weights. 
For this purpose, we must generalize the use of the term "probability" . Instead of a conditional probability, we 
will need what we might term a conditional weight, normalized like a probability so that its integral is unity. 

2 Of course, factorization does not hold as an exact statement. It is only true as an approximation for 
certain important classes of kinematic region. Also, when higher-order corrections are included, we obtain not 
an actually factorizing formula for the cross section but a sum of factorized terms, both for the hard scattering 
and for the parton branching. 
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Figure 3: "Handbag" diagram for DIS. 

The parton-model approximation is obtained from handbag diagrams, as in Fig. ||. Notice 
that we have included a bubble representing the full showering and hadronization of the struck 
quark. Up to a normalization factor, which will not concern us here, the value of the handbag 
is 

Q 2 f d 6 k 



<f>(k,P)D(k), 



(2.4) 



2vr J (2?r) 6 

where $ and D correspond to sums over the correspondingly labeled subgraphs in Fig. ||, 
with an integral over all final states. The parton-model approximation Q applies when k has 
low transverse momentum and virtuality compared with Q and when the struck quark also 
has virtuality much less than Q. Then in the hard scattering, which in this case is a trivial 
factor of unity, we can replace the incoming and outgoing quarks by on-shell momenta 



k ~ (sP + ,0,dr), 
k + q~ (0,g~ T ) = (0,Q 2 /(2xP + ),0 T ). 



(2.5) 



However, this approximation cannot be made in the two nonperturbative blobs $ and D, but 
rather the small components of momentum must be integrated over. This leads to 



xP^ 



dk d kx 
(2vr) 6 



$(xP + ,k-,k T ,P) 



2q-_ 
2vr 



dl + D(l + ,q-,6 T ) 



(2.6) 
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Up to a normalization, the first factor in square brackets corresponds to the usual definition 
of a parton density, f(x), as the matrix element of a light-front operator. The second factor 
is the integral of the fragmentation function D over final states, again up to a normalization; 
this integral is over the discontinuity of a propagator, so that the factor is equal to unity if 
the integral is convergent. 3 

Notice that in the $ and D factors in the parton-model approximation we have replaced 
the large components k + and k~ + q~ of the external momenta by the same approximate 
values xP + and q~ that are given in Eq. ( |2.5| ); it is only the small components that we 
leave unchanged. 4 We have therefore replaced the final states of <£ and D by final states 
with different momenta. For the inclusive structure function F2, this is satisfactory, since 
we are only concerned with the numerical values of the factors in Eq. Q2.6| ), not with the 
detailed structure of the final state. For this case, errors in the approximation at large parton 
virtuality are properly compensated by the usual correctly derived higher-order corrections 
to the hard scattering. 

But the situation is different in a MCEG, as we will now see. In the first step of a MCEG, 
a hard-scattering configuration is generated with a corresponding value of parton momentum 
k. Since the showering has not yet been performed, it is natural and conventional to assign to 
this momentum the naive parton model value, as in Eq. (2.5). The cross section is expressed 



in terms of the parton density f(x), which, as in the first factor of Eq. Q2.6| ), is defined in 
terms of an integral over parton momentum k with the plus component fixed: k + = xP + . 
However the true value is different: 

M 2 + k 2 

k+ =' p¥ + W^' <2 - 7) 

where M is the invariant mass of the outgoing parton. 

Thus there is a mismatch between the true value and the initially generated value of k + 
that is used to define the parton density. This can be a large mismatch, since the virtualities 
and transverse momentum range up to order Q, at least. The parton-model approximation 
replaces the final state in blob $ by a final state with a different momentum. 

Since in a MCEG we resolve the structure of the final state, the value of k has to be reas- 
signed after showering to correspond to the true final state. The use of the word "reassigned" 
in the previous sentence unmasks a conceptual shift that is critical in correctly analyzing the 
algorithm in a MCEG. There are in fact (at least) three very different objects to which we 
can give the name k. Two correspond to the dummy variables of integration in one or other 
of the two formulas ([O]) and (|2.6| ) for the handbag diagram; within a Monte-Carlo imple- 
mentation each of these variables can be treated as a random variable which has a single 



3 There are, in fact, ultra-violet convergence problems in these approximated integrals, associated with the 
need for NLO corrections to the parton model; the divergences can be dealt with either by a cutoff or by 
renormalization But this issue will not affect this part of our discussion. 

4 We have also replaced the parton transverse momentum in the fragmentation by zero. But this is simply 
equivalent to a small Lorentz transformation. 
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definite (Lorentz- vector) value in each event. The third object called k is a storage location 
in the computer program; its value can change during the generation of a single event. Such 
a reassignment of the computer variable changes the correspondence between the computer 
variable and the random variable(s) used to define probabilities. 

Branching probabilities in initial-state showering are conditional probabilities given the 
value of k, and in particular given the value of k + . But the value of k + used in the initial- 
state showering changes during showering, in a way that depends on the shower history 
of the event. The same applies for the value of k~ + q~ in the final-state shower. This 
makes the fundamental basis of a conventional MCEG and of its independent showering 
probabilities hard to quantify accurately enough to allow a transparent derivation of higher- 
order corrections. 

A useful characterization of the difficulty associates it with the momentum-conservation 
delta function at the photon vertex. Let I be the momentum of the outgoing quark. Then 
the standard parton-model approximation replaces 

S^ik + q-l) (2.8) 

by 

5(k + +q+)5(q~ -r)5W(T T ). (2.9) 

These are good approximations to each other only if they are integrated with a sufficiently 
slowly varying function, as is the case for the inclusive cross section. But when we resolve 
the full final state, the two delta functions are not good approximations to each other; they 
are infinitely different. 

Therefore we propose that the fundamental factorization formula for a MCEG should 
involve unintegrated parton correlation functions, and that showering then involve probabil- 
ities conditional on k and k + q which have already been assigned their exact values. While 
the approximated momenta will still make their appearance, it will be in a different and 
controllable fashion. 

By using a factorization formula which does not contain ordinary parton densities we 
have also removed the DGLAP evolution kernels, which appear to form an essential part of 
the showering algorithm. Other quantities will take over their role. 

We will find it convenient not to generate the exact parton momenta at one step, but to 
reformulate in several substeps the first step of the algorithm that we gave at the beginning 
of this section. Thus we can generate the hard scattering as follows: 

l.a Generate the hard scattering with conventional massless external partons, and conven- 
tional parton densities. 5 

l.b Generate, independently, the virtualities and transverse momenta of the partons ac- 
cording to some suitable approximate distributions. 

5 The parton densities should be the most accurate known parton densities, rather than just the approxi- 
mation with LO evolution. 
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l.c Reassign the parton momenta to their exact values according to a defined prescription. 

l.d Reweight (or appropriately veto) the event so that the probability/weight corresponds 
to the correct distribution in terms of parton correlation functions. The values of 
the parton correlation functions will have been tabulated as the result of a separate 
calculation. 

A similar procedure is applied at each step of parton showering. The motivation of this pro- 
cedure is that the factorization formula in terms of unintegrated parton correlation functions 
involves more variables and more complicated objects than standard factorization. So it is 
simpler to base the event generation on ordinary parton densities and then correct the result 
to correspond to the new factorization formula. 

Although the reassignment of parton momenta and the reweighting are superficially sim- 
ilar to those in other proposals, the important point is that we propose a change in the logic. 
Previously the hard scattering and the branchings steps in showering generated approximate 
parton momenta. The reassignment of parton momenta to their exact values occurred after 
showering of these partons. If we then try to use a reweighting method to treat higher- 
order corrections for both the hard scattering and for showering probabilities, there is overlap 
between the parts of the event to which different reweightings are applied. 

We instead do the reassignment of parton momenta before the showering. All the proba- 
bility distributions for the actual parton momenta are distributions as functions of the exact 
momenta. Therefore the branching probabilities for different partons can be taken as strictly 
independent of each other. Moreover the reweighting is not to make the differential cross 
section correspond to the standard NLO calculation of the cross section but rather to make 
it correspond to our new factorization formula in terms of parton correlation functions. Of 
course, the new formula will include the information in the previous calculations, but it 
can potentially include even higher order corrections. Most importantly, the application of 
reweighting to the hard scattering or one branching does not overlap logically with later 
branchings. Thus exactly similar methods can be applied to both the hard scattering and to 
parton branching, so that both can be systematically treated to arbitrary order. 

Our procedure is a generalization of the work by one of us in § for final-state showering. 
The new feature is in a certain sense concerned with the directions of cause and effect. In 
e + e~ -annihilation we generate a hard scattering given the incoming leptons. We shower 
the partons and readjust their kinematics by some chosen prescription to satisfy momentum 
conservation. 

But with an initial-state hadron, this procedure requires us to readjust the kinematics of 
the incoming parton whose value of k + was used to determine the cross section from the parton 
density. The transformation of the kinematics to get correct parton momenta inevitably 
changes the value of k + — see Eq. (2.7). This invalidates the value of the non-perturbative 



parton density used at the initial step of computing the hard-scattering probability. The 
corresponding non-perturbative probability in e + e~-annihilation is the unit probability for a 
parton to shower. 
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3. The model for DIS 



Our method in its present form works for any non-gauge theory. But for simplicity of pre- 
sentation, we work with (ft 3 theory in 6 space-time dimensions. The Lagrangian is 

C = \{d4>? - \m*f - \g (^) 4> 3 + C ct , 

Cct = bz 2 (d^) 2 - Um 2 <f - Ug^ - 8h<p. (3.1) 
2 2 b 

For regularization of ultra-violet divergences we use a space-time dimension n = b — 2e. 
We associate a factor fi e = [ju 2 e 7 /(47r)] e / 2 with the coupling so that MS renormalization 
is implemented by pure pole counterterms. We use MS renormalization for all except the 
renormalization of tadpole graphs, and for those we define the Shtp counterterm || by requiring 
the sum of the tadpoles and their counterterms to vanish, i.e., by the renormalization condition 
(0|(/>|0) = 0. We define a bare parton field as <po = Z l / 2 cj). 

Our model for DIS consists of the exchange of a weakly interacting boson, A, for which 
we add an extra term to the interaction C s = eo^40o/2 and a corresponding counterterm. As 
in HI, we unify the calculation of all cross sections in a weighted cross section 

a[W] = Kj2(P\m\f)(M0)\P) (2nf5^(q + P-p f ) W(f) 
f 

= Kj2W(f) [ d 6 y e 1 ^ (P\j(y)\f)(f\j(0)\P), (3.2) 

where the weighting function W(f) is an arbitrary smooth function of the final state. The 
current j = e^[4> 2 ] = eo^</>o contains a renormalized composite operator H [(f) 2 ], and K is a 
standard leptonic factor. We have chosen to include a factor of the analog of the electromag- 
netic coupling in the operator. This ensures that as regards strong interactions the operator 
is scale invariant; the "electromagnetic" coupling e in our model gets renormalized by strong 
interaction effects unlike the case of the true electromagnetic coupling in QCD. 

Any exclusive component of the cross section can be obtained from a[W] by functional 
differentiation with respect to W, while, by the choice of a suitable form for W(f), jet cross 
sections and in general any kind of binned cross section, can be obtained. The use of weighted 
cross sections is very compatible with the MCEG approach, since the output of a generator is 
a list of events with a distribution corresponding to an approximation to the fully differential 
cross section. An estimate of a[W] is obtained from a sum over the generated events weighted 
by W(f): 

<r e8 t[W] = YJ2 W (M> ( 3 - 3 ) 

i 

where fa is the final-state of event i, and L is the luminosity appropriate for the set of events. 
If the EG gives weighted events, then a factor of the event weight Wi needs to be inserted 
inside the summation. 
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4. Hard-scattering factorization 



The key idea of our method is that all partons should be generated with their exact kinematics. 



So in this section we will obtain a factorization theorem, Eq. ( 4.31 ) below, which expresses 
the general weighted cross section a[W] in terms of parton correlation functions. 

While we will use the new factorization property for generating the partonic content 
of DIS events, we will also need the usual factorization theorem to determine values of the 
inclusive cross section: 

a[l] = K I — C(x/£, Q 2 /^ 2 , <?(/u)) fi) + power-suppressed correction, (4.1) 

Jx ? 

where C is the hard scattering coefficient and / is the ordinary parton density in the target. 
We choose to define the parton density in the MS scheme. 6 The coefficient C is perturbatively 
calculable if is chosen to be of order Q. Notice that this theorem only applies after an 
unweighted sum over final states, so the argument of a is a weight function whose value is 
unity for all states. 

The normalizations are chosen so that the lowest order value of the hard-scattering co- 
efficient, associated with the handbag diagram with neglect of hadronization, reproduces the 
parton model: 

C(x/H, Q 2 /v 2 ,g) = e 2 i8(i - x) + 0(e 2 g 2 ), (4.2) 



so that 



<t[1] =Ke 2 f(x;v) [l + O {g 2 (Q))] 



(4.3) 



4.1 Regions 



P 



ui 
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Figure 4: Decomposition of a graph for DIS corresponding to a generic leading region. 



As usual, the dependence of / on the renormalization scale fi includes the effect of the running coupling 
and mass g(fi) and m(fi), and is governed by a DGLAP equation. 
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The regions that give leading power contributions for a DIS process are obtained by standard 
arguments || and are illustrated in Fig. |j: Momenta in the lower subgraph F are approxi- 
mately in the proton direction, and momenta in the upper subgraph U are grouped around 
one or more different directions. The hard subgraph H has internal momenta that are all off- 
shell by order Q 2 . If momenta could be cleanly separated into different regions, a sum over all 
graphs and over all regions for each graph would give quite directly the factorization formula, 
Eq. fl4.31| ) below, which corresponds to the factorized structure in Fig. ||. Moreover inter- 
preting the blobs in Fig. || as complete amplitudes, not just as sums over Feynman graphs, 
strongly indicates what factorization should mean non-perturbatively in the full theory, even 
if we are not yet able to make a real proof beyond perturbation theory. 

However, there are momentum configurations that interpolate between core parts of the 
different regions and these give well-known logarithmic enhancements, evolution of parton 
densities, etc. So an unambiguous and clean separation between different regions of momen- 
tum does not exist. All but the simplest graphs for the cross section (or structure function) 
have several leading regions. We will analyze them region by region, obtaining an expression 
for each graph T as a sum of terms, one for each region R: 

T = ^ Cr(T) + power-suppressed correction, (4.4) 
R 

with the application of subtractions to compensate double counting between regions. The 
principles of the argument will be exactly those of our previous work [Q] for e + e~ annihilation. 
The changes will be those to accommodate an initial-state hadron and an algorithm for 
backward showering for the initial state. 

We will find that subtractions only need to be applied inside the hard scattering, so that 
summing over graphs and regions will give a factorized form for the weighted cross section. 
Schematically 

a[W] = K R IT ^ R W(f) + power-suppressed correction. (4.5) 

/ graphs T j=l 
regions R 

We use a subscript R to denote the values of quantities that are appropriate to a particular 
region of a particular graph. This will allow us to use the unsubscripted values for the 
corresponding quantities H, and F, etc, summed over the relevant subgraphs, that appear 
in the final factorization formula. Thus Hr is contribution to the subtracted hard scattering 
factor, from the H subgraph in Fig. [|. Similarly Fr and Unj correspond to the target and 
final-state "jet" factors. We have made explicit the sum/integral over final states, so that we 
can insert a non-factorized weight W{f) to analyze the final-state structure. We will see that 
the sum over graphs can be performed separately for Hr, Fr and Urj, and this will lead to 
the actual factorization formula, Eq. ( [4.31 ) below. 
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4.2 Factorized approximation for one region 

We now construct a factorized approximation for a general leading region R of a general 
graph r for DIS. The approximation must be compatible with the structure of a MCEG, 
as described in Sec. |], which means that it must have a factorized form and must preserve 
exact parton momenta for the final state. The accuracy of the approximation will degrade 
as momenta move away from the core of the region, and the full factorization theorem will 
contain subtractions in the definition of the hard scattering coefficient that eliminate double 
counting when we combine the contributions of different regions. Then the final answer will 
be uniformly accurate in all the different regions and in the interpolating regions. 
It is useful to define a systematic notation for the momenta in Fig. 

• N u is the number of the final-state partons connecting the hard scattering and the upper 
"jet" subgraphs in U(q, k, m). 

• rij is the number of final-state particles arising from parton j in subgraph Uj. 

• The momenta of these particles are pj t i, pjp, . . . , Pj, nj - 

• p., j = 1, . . . ,N U , with an underline, denotes the collection of the final-state particle 
momenta in the jet subgraph Uj, i.e., p. = (j>j ; i, Pj.2, ■ ■ ■ > Pj.n,)- 

• lj = Xyi=iPi,i denotes the momentum of parton j, and Mj denotes its invariant mass. 

• l u = (Ji, I2, ■ ■ ■ , ljsr u ) denotes the collection of the final-state parton momenta. 

• M_ = (Mi, M2, . . . , Mj\r u ) denotes the collection of the final-state parton invariant 
masses. 

• P = (Pi j Pn, ■ ■ ■ ) P M ) denotes the collection of all the final-state particles in the upper 
part of the graph, U(q,k,m). 

• Similarly we use index j = for the target subgraph F, i.e. 

— uq is the number of final-state particles in the target subgraph F. 

— p Q denotes the collection of the final-state particle momenta in F. 

We introduce the following notations for the final-state phase space: 



dL(p ,q + k; mph ) = Jin ^ I (^)^ {6) (E^ 

d 5 Pj, 

1 (2vr)52^ + m ph 



k), 



dLip-lj-m^) = [] (2vr) 6 ^ 6 )(E^-/,), 



V d 5 L 



dL(l_ u ,q + k;M) = TT i (2^ (6) (E jlj ~ Q ~ *0- (4-6) 

fi (27r)52 A //l + M, 2 
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Here m pn is the physical mass of the true final-state particles. As pointed out in ||, the 
integration over the final-state phase space is factorizable: 

dL( gu , q + k;m ph ) = nflT^lk (^^Ei.^,* " 9 - *) 
j=i i=i ^'■ t 

n/^/ ^ ,) ( 2^<«>(E ii3 - 9 -*) 

IX 

n 



-AT, 

X 



,i=i 



ft Mm?\ /■ ft / \ 

n y -^r i y ^Ou^+^Mjn^fe^ 7 ^) • ^ 

With this notation, we now write a decomposition for a graph T and a particular region R of 
the form in Fig. ^ 

r[VF] = / (Se dFR^k.P) W(f) (4.8) 

with 

dF R (p ,k,p') = dL^P-k-m^ F R (p Q ,p) . (4.10) 

Note that at this point we have not changed or approximated the graph in any way; we have 
simply written the graph in a form suitable for analyzing it in region R. 

In 0, integrated jet factors were defined corresponding to sums of each of Urj over graphs 
and final states. These serve as normalization factors in the MC algorithm for showering. 
Because of the different organization of the algorithm in this paper, our definitions will be 
slightly different. First we define an integrated jet factor as 

J/(/ 2 ;/i) = / d?y <O|0(2MO)|O) e*». (4.11) 



For consistency with our naming for the target-related factor, we can also call this a "vacuum 
parton correlation function" or a "final-state parton correlation function" . 

Similarly, for the target factor, we define a parton correlation function (pCf) §>(k,P) as 
the full non-perturbative quantity that corresponds to integrating F over all final states and 
then summing over graphs: 



*(k,P)= Yl dL(p ,P-k;m ph ) Fr^P^ 

graphs 
final states 

= J d e y (i>|0(j/)0(O)|P) e~' lk -y. (4.12) 
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The pCf is necessarily nonnegative: 



Y ^ 

= J2(^fS(k + Py-P) \(P\m\Y)\ 2 



Y 



> o, 



(4.13) 



as is the integrated jet factor Jj. 

We could now formulate factorization with an off-shell hard scattering. But on-shell 
and massless amplitudes are easier to compute. So we will now define an approximation 
appropriate to region R in terms of a hard scattering defined with massless external partons. 
We will denote the approximation by an operation Tr: 

T R T[W] = ]^J J dL(p ,P-k,m ph ) F R (p Q ,P) 



W(f). (4.14) 



x / dL(l_ u ; k + q; M = 0) H R (l_ u , fc, m -► 0) 



n 



In this formula we have used a projection of the true parton momenta k, h, . . . , Zjv„ onto 
massless momenta k, h, . . . , Zjv u) with zero transverse momentum for k; we will define such 



a projection in Sec. 4.2. We have defined the hard scattering to be computed with massless 



external momenta and with zero mass parameter for the internal lines. However, the target 
factor Fr and the final-state jet factors Ur. are computed with the exact parton momentum 
k and l u , not the approximation. 

As in Q, there is a cut-off function Q(Mj / /j/j) that restricts the final-state parton masses 
to be below some factorization scale 7 fij, which should be of order Q. An ordinary ^-function 
6(M?//Xj) = 9(1 — Mj/fij) would be suitable, although a smooth cut-off might be better 
for implementation in MCEGs. After we construct a factorization formula, this cutoff func- 
tion will also appear in the subtractions in the hard scattering. Since the cutoff function 
does not appear in the definition of the cross section itself, there is a kind of generalized 
renormalization-group invariance: Changes in the value of \xj and in the functional form of 
will leave the factorized approximation to the cross section unchanged, up to power-law 
corrections, when all orders of perturbation theory are included. 

A similar cut off function could be applied to the initial state, but the coupling of the 
kinematics of initial and final state showers (through energy-momentum conservation) com- 



7 Conventionally, the term "factorization scale" applies to a quantity analogous to fij but for the treatment 
of initial-state parton radiation in the conventional factorization formalism. But the same principles apply in 
fragmentation functions and the like in the final state, so we will use the same name. In ||, the symbols /if 
and ur were used to denote the quantities that are {ij and [i in this paper. 
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bined with our use of exact parton kinematics from the beginning for the struck parton's 
momentum k, removes the need to do so. 

Changing integration variables to the exact parton momenta gives 



T R T[W]=l[J JdL(p o] P-k,m ph )F R (p ,k,P)H R (i u ,k,m^O) A Nu (l_ u ,k) 



n 



dL(p • Z i; m ph )G(M 2 / n 2 j)U Rj (p 



W(f), 



(4.15) 



with k = ^2jlj — Q- The quantity An u is the Jacobian of the transformation between the 
exact momenta and the massless approximate momenta — see Eq. ( [4.22 ) below for its value. 



In a region where the final-state lines of the subgraphs U R j do in fact form jets and where 
\k 2 \,k^ <C Q 2 , T R F provides a good approximation to the original graph T, up to power-law 
corrections (and the Jacobian approaches unity). 

Now when we treated e + e~ annihilation in Ref. Q we found a convenient and simple 
factorized approximation in terms of the massless momenta. It allowed the MC generation of 
massless parton configurations for the hard scattering, with probabilities independent of the 
showering. But in DIS, the coupling of the kinematics of final and initial-state partons, as 
given in Sec. ||, prevents us from having a correspondingly simple decoupling of showering from 
the generation of massless partons for the hard scattering. So it is not so clearly advantageous 
which set of variables to use here. It follows that there is also a choice of which formula for 
T R T is defined to contain the Jacobian Ajv u . Any change here is compensated for in physical 
cross sections by the subtractions in higher-order corrections to the hard scattering, so this 
is primarily a matter of convenience; the Jacobian is unity in the collinear region. Our choice 
probably simplifies the logic of deriving subtractions for hard scattering, since they most 



naturally involve massless approximations at certain points. See Sec. 4.7 below for further 
details. 



4.3 Definition of projection onto massless momenta and its Jacobian 

In this section we give one possible definition of the projection onto massless parton momenta: 
k = PN u (k), l u = Pn u (L u )- We require it to obey 

k~ = k T = 0; (4.16) 
/] = 0, (j = l,...,N u ). (4.17) 

We also require energy-momentum conservation, so that Y2j h = 1 + ^- This gives the 
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following constraints; 



i * 2xP+ ' 

3 3 

-hr + ^ l iT = ^2 ^ T = Qt = 0; 

3 3 

-k + + J2 l t = + Yj*t = i + = ~ xP+ - ( 4 - 18 ) 

3 3 

In addition, this projection should have a smooth limit when some or any of the external 
momenta go to zero. Any projection would be a valid choice if it satisfies the above constraints 
and Pn u [Pn u {1, k)\ = Pn u (^ k). It is also desirable that the formulas be invariant under 
boosts in the z direction. 
Our chosen definition is: 

(k T = 0, 

VjT = l jT ~ k T /N u , 

(k- = o, 

\lj = IJ -k~/N u , 

l) T _ (Ijt ~ k T /N u ) 2 



3 217 2(17 - k-/N u ) 



U 3 U2{l--k-/N u ) 



(4.19) 



Given the values of A; , kr and of the Mjs, which parameterize the deviation of the true 
parton momenta from the projected momenta, we find the inverse transformation: 

l jT = l jT + k T /N u , lj = lj + k~/N u , 

I 2 4- M 2 

n = ^^, fc+ = (4-20) 

3 3 

The kinematics of DIS requires that k 2 < 0, so either a rejection or a good selection method 
is needed to eliminate positive k 2 events. 

Notice that, although DIS kinematics require that qx = 0, our choice of projection applies 
to qx / without any modification. 

We need the Jacobian of the transformation between the exact and the approximated 
parton momenta. We obtain it by expressing dL(l) in terms of — and T components of 
momentum, and then using the fact that at fixed k~ and kx, the transformation from lj and 
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IjT is a simple shift. Thus 



J "« - r dh d% T 1 I sr^ j 



(2tt) 5 2/J- 




dlj d IjT 


1 


(2n) 5 2lJ 


2 W 


dlj d IjT 


1 


(2vr) 5 2f- 





AT U -1 

/ AT U I- ' 

dk + dL(l u ;q + k,0) \ Y[j^)- (4-21) 




Hence the Jacobian is 



A (1 h\ - dk + dL(l u ;q + k,0) _ yr l j _ yr l j , , 

A ^ k) = dk +d L(l u ; q + k;M) " H f " H I; - (422) 

In the collinear limit for the initial-state, where k~~ and &t approach zero, the Jacobian 
approaches unity. 

4.4 Construction and proof of factorization 

The proof of factorization proceeds as in P], with very few changes. We first define the 
contribution Cr(F) associated with a leading region R of a graph T, using the operation Tr 
and some subtractions. Then we show that with this definition, the sum over R of Cr(T) 
gives a good approximation to the graph everywhere, i.e., the remainder in Eq. ( [4,4]) is power 
suppressed. Then we sum over graphs for the cross section. Using the structure of TrT and 
of Cr(T) we show that this gives a factorization. For the most part it will be sufficient to 
indicate the differences from [||. 

To allow factorization to be a practical tool for calculations, we need the hard-scattering 
factor H to be perturbatively calculable, and we need to construct an algorithm for showering 
with perturbatively calculable evolution. Perturbative calculability of H follows just as 
since we have the same subtractive structure; if we set the renormalization and factorization 
scales [i and \ij to be of order Q, then the integrals for H are dominated by hard momenta on 
a scale Q. We will construct a suitable showering algorithm in Sec. ||. Before we do this, we 
will, in Sec. |j| analyze the properties of the parton correlation functions; this analysis re-uses 
the methods for factorization of the cross section. 

Term for a given region As in [g], each leading region is labeled by a particular pinch- 
singular surface (PSS) § in a massless theory, and, as illustrated in Fig. [|| a leading region 
R of a graph T is characterized by a decomposition into subgraphs: a hard subgraph, one 
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initial-state "jet subgraph", and one or more final state "jet subgraphs". We defined an 
approximation T R appropriate to the region in Eq. ( |Og ); its action can be symbolized as 
converting 

T = F R xH R xU R1 xU R2 x---x U RNu (4.23) 

to 

T R T = F R x {TrHr) xU Rl xU m x---x U RNu x A Nu . (4.24) 

For the representation of the graph as a sum over contributions from regions, we will define 
the contribution of region R by applying T R to the graph, but only after subtracting the 
contributions from smaller regions, to avoid double counting. As in |^], we slightly modify 
this, to include a "veto factor" V R : 

c R (r) = v R T R I r - Y, C R' T ) • ( 4 - 25 ) 

V R'<R ) 

These equations give a recursive definition of the contribution C R (T) for region R, which 
starts from the smallest regions. For a minimal region, Rm'm i-6. ; one with no subregions, the 
definition is simply C Rinin = T Rinin F. 

Veto factor The veto factor arises because @ the integrals in the cross section include 
regions in the hard scattering where partons are collinear. Now the subtractions in the 
hard scattering imply that the corresponding contributions to the cross section are power 
suppressed. Thus in infra-red- and collinear-safe cross sections, the veto factor is unnecessary. 
But the contributions are not completely well behaved, so we define a veto factor, to eliminate 
small regions around the singular surfaces for the massless hard scattering, as follows: 



Vr = H \v(M 2 s /m 2 ) V{-~kl/m 2 )] . (4.26) 



s 

The product is over all subsets S of the (massless) final-state parton lines of the hard scat- 
tering. In this formula Ms is defined to be the invariant mass of these lines, i.e., 

while k s is defined to be the momentum transfer with respect to the incoming parton line: 




_ (4.28) 

The elementary veto factor is defined as 

V(M 2 /m 2 ) = 0(M - m). (4.29) 
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Thus the veto factor Vr for region R eliminates a small neighborhood of the initial-state and 
final-state collinear singularities associated with the massless partons in the hard scattering. 
The exact definition is irrelevant. 

We will show an example of the initial-state veto factor in Sec. [/]. 

Factorization Now that appropriate definitions have been made, the proof of factorization 
given in Q applies to the case of DIS. When we construct Cr(T), the subtractions are for 
smaller regions. As in Q , a smaller region implies larger collinear subgraphs and therefore a 
strictly smaller hard subgraph. Therefore the subtractions are applied only inside the hard 
subgraph: 



Cr(T) = F R x V R T R (h r - Y, Cr'Hr] x[/ ffl x • • • x U RNu x A Nu . 
V R'<R J 



(4.30) 



This contrasts with the case of gauge theories, where the presence of leading effects due to 
soft gluons gives a more complicated situation, which needs disentangling with the aid of 
Ward identities. 

The factorization formula is now obtained by summing over all regions and graphs. The 
correspondence between regions and subgraphs converts the sum to independent sums over 
graphs for Fr, Hr, N u and the Urjs. Hence 



r Nu r i 

a[W] =KJ2 Hn u A Nu F(p Q ,P) n [J(PJ W(f), (4.31) 

3=1 



up to power-suppressed corrections, as usual. The sum and integral represent a sum and 
integral over final states and then a sum over decompositions of the partonic structure and 
over the assignment of final state particles to the F and J factors. Explicitly, the factors are 
as follows: 



1. The hard factor Hn u is obtained by summing the hard scattering factor in Eq. ( 4.30 ) 
over all graphs for Hr: 

H Nu (kJu;,Q,Vj,v,g(v)) = E VT [ H R- E Cr'Hr) > ( 4 - 32 ) 

graphs \ R'<R / 

with a given number of final-state partons. The relevant graphs are for an amputated 
amplitude times complex conjugate amplitude with N u final-state parton lines, one in- 
coming parton line, and one virtual photon. The unsubscripted T denotes the operation 
of projecting the external partons onto massless momenta and setting to zero the mass 
m on internal lines, while V without a subscript is the corresponding veto factor. The 
proof in shows that H has no collinear divergences. 
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2. The initial-state factor F is the sum over graphs of Fr for a given final-state. Thus we 
have an operator formula 

2 



F(p ,P;v)= (P\<P(0)\P 

The integral of F over all final states with a given value of k (k 
the pCf 



/ dL(p Q , P-k, m ph ) F(p , P) = *(fc, P). 



(4.33) 
EiPo.i) gives 

(4.34) 



3. Similarly we have a final-state factor for each final-state parton, which we define exactly 
as in 101: 



It obeys a sum rule 



E 



J(p,;/i)= (O|0(O)b 



dLip.^l^m^) J(p.) 



(4.35) 
(4.36) 



We now reorganize the factors in Eq. (4.31) to be suitable for a MCEG. We assume 
that we have generated values of the usual x and Q from standard factorization ( f4.1[ ), these 
variables determining the lepton kinematics. The problem is now to find probability densities 
for the internal variables so that the generated events give correct weighted exclusive cross 
sections. So we first normalize the factorization formula for a[W] by dividing by the inclusive 
cross section cr[l]. Then we choose as independent variables the massless variables, k and 
lj, together with the variables k~ , kx, and Mj that characterize the deviation from exact 
collinearity. We write a factor that can be used to generate a hard scattering in terms of 
the massless projected momenta independently of the showering. Then we write a factor 
that gives an approximately normalized 8 probability density for k~ , hr, and factors that give 
normalized probability densities in Mj, and for the final states resulting from each final-state 
parton. Finally we insert a factor that expresses the deviation of the true formula from an 
exactly factored form. We find 

dk+ dk- d 4 k T 



a[W] 



oo 

E 

N u =l 



(27T)' 



da Nu {Q,fij,fi)_ £ + + _ $((fc+,fc ,k T ),P;n) ^ 



a[l] 



f(k /P + ;v) 



X 



J'=l 



n 

=1 

E 



"0 



»dMf Ji(Mj;n)Q(Mj/ri) 
dL (P , P-k; m ph ) F(p , P; //) 



$(k,P;n) 



H(k + ,k-,k T ),P;p) 



$(k,P;v) 
W(f). 



N u 

n 

i=i 



E 



dD{p.;lj;fi) 



(4.37) 



In a sense to be discussed below. 
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Here / is the ordinary renormalized parton density, discussed in Sec. 5.2 below. Notice that 
we have defined by an operator matrix element in the exact theory, in which the particle 
mass is nonzero. Thus the integral over k is free of the collinear divergences that arise in 
conventional calculations in a massless theory. 9 In addition, k is the exact quark momentum, 
so that the upper limits on its components are set by the kinematics of the process, and there 
are no ultra-violet problems. Thus the integral is convergent. 

In the first line we use what we will call the differential hard-scattering cross section 
defined by 

da Nu =K dL(i u ,q + k,0)H Nu (Q, t ij, t i)Z{^j/fi,m/n) N -. (4.38) 

As in Q , we have inserted a factor Z for each outgoing parton. This provides a kind of analog 
to the propagator residue factors in the LSZ reduction formula, and is defined by 



00 dM} 
» dM 2 



Zijij/^ m/fj) = I e(MjV/4) J/(M?; /i) (4.39) 



~2^^J dL (P r l j)®(Mj/»j)J(p r ri- (4.40) 

Observe that j|] since the integral of J/ over all Mj diverges at large Mj, we have to use a 
cutoff, the factorization scale fj,j, to define a suitable finite quantity. We have implemented 
the cutoff by the function 0(M? / (£j) rather than by the upper limit of integration, to allow 
for the possibility of a smooth cutoff. 

But there appears to be no unambiguously appropriate analog of the factor Z for the 



incoming parton. Instead, on the first line of Eq. (4.37), we have multiplied da by the 
ordinary parton density f(k + /P + ; fi) taken at the approximated parton momentum. This is 
like the parton density times hard-scattering in the conventional formalism, and we can use 
the numerical values of daf for generating the kinematic variables k + and lj. However, as 
we have already argued, the value of k + does not agree with the correct value of true parton 
momentum. Since the correct value cannot be known algorithmically until the full parton 
kinematics are known, we insert a correction factor on the last line of Eq. ( 4.371 ) to adjust 



the full cross section to its correct value and we will implement the correction factor as a 
reweighting. 

To a first approximation, the incoming parton for the hard scattering has low transverse 
momentum and virtuality, and then k + is close to k + . So daf/a[l] is approximately a 
normalized distribution in k + and the outgoing parton momenta. Hence the reweighting 
corrections are mild. 

The remaining factor on the first line, $//, is therefore also approximately a normalized 
probability density in k~ and kx, with upper limits determined kinematically. The density /, 
conventionally called an "integrated parton density", is therefore approximately an integral 



9 We only make a massless approximation in the fully subtracted hard-scattering coefficients, where the 
approximation is valid. 
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dD(p^) = - 3 in2 _- ] (4.41) 



over However, because of complications from ultra-violet behavior these statements are not 
exact; the detailed relations between $ and / will be explored in Sec. ||. Again, reweighting 
corrections to compensate for slightly inaccurate normalization are mild. 

In the second line, the distribution in outgoing parton mass is given by the Ji/Z factors. 
Given our cutoff at Mj = /Uj, this is a correctly normalized conditional probability density. 

The normalized differential distribution of final-states in an out-going parton given the 
parton's exact momentum is the quantity dD defined in Eq. (6.42) of Pj. In our current 
notation it is 

J(p.;n) dL(p.;lj) 

The corresponding normalized initial-state factor is dL(p Q )F/$. 

The deviation from the factorized form of the probabilities is given by the ratio of the pCf 
at the correct incoming parton momentum to the value at the approximated (but off-shell) 
momentum, i.e., by the factor 3>(fc, P)/$((fc + , k~, h?), P)- This factor approaches unity in 
the collinear limit. In a MCEG it can be implemented by either weighting generated events 
by the factor or using a suitable veto strategy. Because the weighting factor is well behaved 
in the collinear limit, the veto method is reasonable. The important point for the logic is 
that this reweighting factor is used only within the generation of the hard scattering, which 
in our method includes the generation of the exact off-shell momenta of the external partons 
of the hard scattering. 

4.5 Choices of \i and fij 

The renormalization scale \x and the factorization scale \ij are artifacts of our method of 
treatment of the theory, so that the cross section does not depend on them. As usual, 
we exploit the dependence on these variables of factors in factorization theorems to enable 
perturbative calculations to be done without large logarithms. Scale dependence cancels in 
predictions of cross sections up to the effect of uncalculated higher order corrections. 

First, the inclusive cross section is calculated from the standard factorization formula Eq. 
( |4.l[) with /i of order Q; then the coefficient function is a weak-coupling expansion in g(fi) 
without large logarithms. The parton density at scale /i ~ Q is related to the parton density at 
a fixed scale Qo by use of the ordinary DGLAP evolution. The evolution is perturbative. Our 
approach is to define the parton density by MS renormalization, so only the renormalization 
scale /i appears here; there is no separate factorization scale. 

Next in the generation of an event is the use ofEq. (p7|) to give the kinematic distribution 



of the hard scattering. Here we also set // and /ij of order Q so that the hard scattering is 
perturbative. For consistency we must use the same values everywhere in Eq. ( [4.37| ). The 
pCf is an operator matrix element, so it depends on jj, but not jjLj. Its calculation and 
scale-dependence will be considered in Sec. || 

The treatment of the distribution of parton mass Ji/Z, on the second line of Eq. ( [4.37 ), 



is exactly as in M . Since \i and /i j are the same as in the rest of the formula, they are of order 
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Q. The calculation of Z(hj/ f/,,m/ (a) is then perturbative, in powers of g(fi), with the mass 

f?; /x) will have logarithms of Mj t 



dependence being negligible. But J/(M 2 ;/x) will have logarithms of Mj/fi, which get large 10 



when Mj < Q. Its evolution is governed by the anomalous dimension of the parton field and 
therefore by the same RGE as the pCf, Eq. ( |5.2| ) below. Thus we use the renormalization 
group to express it in terms of Jj at another scale: 



J r (M|;/i) = exp 



V dv 

—i(g(v)) 

i V 



Ji{M]-v!). (4.42) 



When Mj is in a perturbative region, we set // of order Mj, so that //(M 2 ;//) can be 
computed perturbatively in powers of <?(//). When Mj is too small, we set // to a fixed scale 
Qo and use a non-perturbative model (fit to data). The model is constrained by a sum rule 
that Ji integrates to the perturbatively calculable quantity Z. 

In the third line of Eq. flO?] ), we have factors that give the distribution of hadronic states 



given the exact parton momenta. Each factor is a phase-space differential times a ratio of two 
matrix elements. The matrix elements differ only by whether or not they are integrated over 
final states, so the numerator and denominator have the same anomalous dimension. Thus 
we can replace fi separately in each factor by a scale appropriate to each parton: of order Mj 
for a final-state parton or of order y^A; 2 ! for the initial-state parton. (A fixed scale would be 
more appropriate in the non-perturbative region of momenta.) Then we exploit factorization 
of the pCf to further the calculation. The change in the value of the scale /i, potentially 
very different in different factors, is the key feature that enables the MCEG approach to give 
systematically reliable estimates of cross sections. 

The final factor &(k, P) /&((k + , k~ , k?) , P) is again a ratio of quantities with the same 
anomalous dimension, so that it is scale-independent. The values needed are obtained from 
the calculation of the pCf. 

4.6 Collinear safety, renormalization group equations 

For the factorization formula to have predictive power we need the following: 

• The remainder, i.e., the difference between the factorized formula for the cross section 
and the true cross section, is power suppressed. 

• Infrared and collinear safety of the subtracted hard scattering factors H. 

• Infrared safety of the integrated final-state jet factors. 

• Renormalization-group-like equations for the scale dependence of the jet factors. 

• Theorems for the parton correlation functions. 

The power suppression of the error and the collinear and infra-red safety follow exactly as in 
[0], and the evolution equations for the final-state jet factors are the same. Then as explained 



3 Note that as a simple operator matrix element Ji has no fij dependence. 
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in the previous section, Sec. 4.5, perturbatively based calculations can be made by combining 
suitable choices of \i and [ij in each factor with the evolution equations to relate the factors 
to their values with a common choice of scale. 



4.7 Factorization schemes 

We have a choice in the algorithm for generating parton configurations for the hard scattering: 



(a) We can follow the structure of Eq. ( 4.37 ): (i) Generate a conventional massless hard 
scattering, (ii) Generate values for hr , kx, and the Mj's. (iii) Compute the exact parton 
momenta k and l u . (iv) Reweight by the ratio of the correct pCf to the approximated 
pCf. 

Or: 

(b) We can generate the hard scattering according to a formula in terms of the parton 
correlation functions alone with exact parton momenta. 

In both cases initial-state and final-state showering are performed conditional on the exact 
parton momenta, so the change in algorithm affects only the algorithm for the hard scattering, 
not the rest of the MCEG or the probabilities. For the second case, we can exhibit the 
necessary probability densities by rewriting Eq. ( 4.371 ) in a form closer to that of Eq. ( |4.3l| ): 



a[W] 



oo N u 

En 

JV„=lj=l 



E 



"0 



x W(f), 



KH Nu <S>(k,P)A Nu (l u ,k) 



(2tt)< 

dL(p Q ,P-k;m ph ) F(p Q ,P) 



N u 

n 

3=1 



a[l] 

Yj / dD (Pj' l j^) 



(4.43) 



with k = Y^j lj — Q- 

Which of these algorithms will be better in practice is not yet clear. That will depend 
on practical experience, at least. 

But clearly, if we use the second method, it would appear more natural to remove the 



Jacobian factor Atv„ from Eq. ( 4.43 ). If we did this, then consistency between the two event- 
generation methods requires us also to insert a factor of 1/ A^ u in the formula ( 4.37|) used to 
specify the probabilities in the first scheme. We would also need to remove the factor Ajy u 
from Eq. Q4.31 ). Consistency with the derivation of these equations requires a corresponding 



change in the definition of Tr: we would insert a factor 1/Ajv u in Eq. (4.14) and remove a 
factor Atv u from Eq. ( [4.15| ). 

Since Tr changes, we would find that changes in the higher-order terms in the hard- 
scattering coefficients occur, as is forced by the subtractions in the recursive definition ( p~25| ) 
of the contribution Cr of an individual region. This gives us two schemes, each more natural to 
one of the two event generation algorithms. As usual, changes in the scheme are compensated 
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in higher-order corrections to the hard scattering, so that the physical cross sections are 
invariant up to uncalculated corrections of yet higher order. 



5. Parton correlation functions (pCf's) 

To build an EG for DIS we will need the values of the pCf <£(/c,P), to use in the hard- 
scattering factorization formula. At large k 2 , this will involve a new factorization formula 
that expresses it in terms of ordinary parton densities (defined, say, in the MS scheme). This 



is exactly analogous to the standard factorization property for the cross section, Eq. fl4.1|) . 
The showering of the parton, and hence the analysis of the final states inside will need a 
second factorization formula analogous to Eq. (E~37l ) or (pU3p for the cross section. 



The second factorization formula for the pCf concerns a quantity $[W] whose relation to 
the pCf $ is like that of the weighted cross section a[W] to the unweighted DIS cross section 
cr[l]. To be precise, we define 

$ m (k,P;fi)=J2 I dL(p Q ,P-k;m ph ) F(p ,P; f i)W(f), (5.1) 

n J 



where F is the same as in Eq. ( |4.33| ), Of course, the previously defined unweighted pCf, 
Eq. ([04]), is the weighted pCf with unit weight for all states: <£(&;, P) = <f>[i](£;, P). 



This factorization formula for the weighted pCf implies an integral equation for the 
unweighted pCf in terms of itself, so that it is not an easy equation to solve for $. However 
the first factorization property for the pCf gives the unweighted $ in terms of ordinary 
parton densities, and it therefore gives an effective method for calculating its functional 
form. Given this information, the second factorization is perfectly adapted to analyzing the 
partonic content of the final states, basis for the showering of an initial-state parton. 

In addition we need renormalization group equations so that the renormalization scale /i 2 can 
be chosen of order \k 2 \ to exploit the factorization theorems perturbatively at large \k 2 \. 

The uses of the factorization theorems at large | k 2 \ are perturbative, aided by renormalization- 
group transformations. But at low \k 2 \ in an asymptotically free theory, we must resort to 
non-perturbative methods, modeling and simple data fitting, just as in current MCEGs. 

5.1 Renormalization group 

Let us first observe that the pCf (weighted or unweighted) arises from matrix elements of two 
renormalized fields, so it has a renormalization- group equation 

$(fc,p; A t) = -7( g (/i))$, (5.2) 



with the conventions of Q. As usual, the total derivative applies to both the direct [i- 
dependence of the Green function and the /x-dependence of the running coupling and mass. 
This formula enables us to transform the pCf from /j, ~ Q as used in hard scattering factoriza- 
tion to the value ~ | k 2 | as used when calculating the value of 3> from the first factorization 
or when applying the second factorization for showering. 
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5.2 Standard parton densities 

The definition of the standard MS-pdf in (ft 3 theory is given in & [l(J. Here we review it and 
some of the properties of the pdf. First there is the unrenormalized pdf: 

f (x)=xP + j dy-e- lxP+ y~ <P|<£ (0,y- T )<M0)|P) c . (5.3) 

The subscript "c" on the matrix element means that we take only the component in which 
the fields are connected to the target particle: the vacuum expectation value is subtracted 
out. Implicitly, /o has dependence on the parameters of the theory, <?(//), m(/j>) and /i, as well 
as on the UV regulator parameter, which we take to be the dimension n of space-time. Notice 
that we use bare fields in this formula, which implies that there is an ultra-violet divergence 
caused by that in Z2. But, if in an attempt to obtain a UV-finite renormalized parton density, 
we simply replaced the bare fields by renormalized fields, that would not work because there 
is a further UV divergence caused by the unrestricted integral over the transverse momentum 
and virtuality of the external parton. So the MS pdf is defined to have all the UV divergences 
renormalized. It has the form 11 

f x dvj 

f{x;n) = lim / —K(r],g,n)fo(r)x) 

00 

K( V , g, n) = 5( V -1)+Y,(n- ^~ N K N (v, 9), (5.4) 

N=l 

where K(rj, g, n) is a series of pure poles at the physical space-time dimension n = 6. 

Notice that in the absence of these divergences, the bare parton density would be an 
integral over the pCf with a factor Z2: 

/„(*) = Z 2 J T $(k,P). (5.5) 

Since we can analyze issues of UV behavior perturbatively, this equation suggests that a 
suitably improved analysis would relate / to an integral over $. Hence, the factor <l?/f in Eq. 



(4.37) is approximately a probability distribution in k~ and kf conditional on the value of 
k + . Because of the existence of higher-order corrections to the UV coefficient functions in 
relations between $ and /, we must expect such statements not to be exact, however. We 
will give a detailed analysis in Sec. |5.4| , and merely observe here that, in a leading logarithm 
approximation, we can approximate renormalization by a cut-off. Thus we can approximately 
write the usual pdf at a large scale M as: 

f(x; g(M),m, \i = M) = f %^ Hh P; M) + 0(g{Mf). (5.6) 



n The necessary dependence of the bare parton density on the the regulator parameter n, both directly and 
through the n-dependence of the bare coupling and mass, has not been indicated explicitly in this formula. 
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It is only in this approximate sense that the usual designation of / as an "integrated pdf" is 
justified. 

Observe that the scale dependence of the parton density arises from two sources: the 
explicit cutoff on the k integral, and the anomalous dimension of the renormalized parton 
field in the operator definition of <I>. 

5.3 Kinematics of pCf 

In non-gauge theories, by Lorentz invariance, the parton correlation function &(k,P;[i) de- 
pends on k only through the scalar variables k 2 and k ■ P. Kinematically <1> has the same kind 
of variables as a DIS structure function, so it is convenient to use variables analogous to the 
standard x and Q 2 of DIS: 



The kinematics in the factorization formula (4.37) for DIS correspond to the first graph in 



Fig. ||, so k is space- like. Hence, even though the pCf is defined for all values of R 2 , it is 
only used when R 2 is positive. Just as in DIS, the positivity of the squared mass of the final 
state implies that < £ < 1, when R 3> m. Since the kinematics of the pCf are completely 
analogous to those of DIS structure functions, so are the factorization properties, essentially 
obtained by the replacement of currents by parton fields. 

In gauge theories like QCD, the gauge-invariance of the theory requires us to insert 
Wilson lines in the operator definition of pCf's. Appropriate direction(s) for the Wilson 
lines are those compatible with a proof of factorization, a proof that would be much more 
complicated than in our (ft 3 theory. If is the direction of the Wilson line, then there are 
two extra scalar variables: (P ■ n) 2 /n 2 and k • n/P ■ n. (Since $ is invariant under scaling 
of n only variables invariant under this scaling can be used.) The extra two variables will 
increase the computational cost: for example in a much larger interpolation table to be used 
in the numerical evaluation of the pCf. In addition there will be more complications in the 
factorization theorems for the pCf. Note that it should be possible to obtain the dependence 
on one of the extra variables, probably (P -n) 2 /n 2 , by a suitable generalization of the Collins- 
Soper equation [11]. 



5.4 Relations between ordinary parton density and parton correlation function 

There are two kinds of relation between the ordinary pdf and the pCf The first gives a 
result for <I> when the parton virtuality is large. The second result gives an integral of <1> over 
parton momenta k. 

5.4.1 PCf at large \k 2 \ 

The result for large \k 2 \ arises because, as we have already stated, the pCf is like a structure 
function but with different operators. When \k 2 \ is much larger then some typical nonper- 
turbative scale, it follows that there is a factorization theorem for the unweighted pCf of the 
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same form as Eqs. (4.1) for the DIS cross section or structure function: 



f 1 d£ 

<fr(k,P;fi)= / — C$(£/£, —k 2 /^ 2 , g(n)) /(£; ^) + power-suppressed correction. (5.8) 

k ? 

Calculations of the coefficient C$ are completely analogous to calculations of the structure 
function F 2 , except for details due to differences in the graphs. For example the lowest-order 
non-trivial graph for a pCf of a parton target Fig. || is 

Hk,P^) LO = 2kS((P - kf - m 2 ) { J^ 2)2 = 2^S(C - l) (fc /^ 2)2 - (5-9) 

Notice that unlike the case of /, which is integrated over k 2 , the graph without any interactions 
does not contribute here since we will make \k 2 \ large. The lowest order term in the coefficient 
C$ is therefore 

C*(C/£, R 2 ; fi, g) = ^f(S(C -0 + 0{g A ). (5.10) 
Hence the lowest order estimate of the unweighted pCf is 

Hk, P; fi) = ^0^f((; /x) + O x logarithms) , (5.11) 

where, as already mentioned, £ = k 2 /(2k ■ P). Note, that this 
formula is only to be applied at large \k 2 \. At low k 2 we cannot 
neglect the parton mass, and therefore the apparent singularity 
at k 2 = in Eq. ( 5.1 1| ) is not physical. 



ky 









In the higher terms in Eq. ( 5.11[ ), there are, as usual in 



Figure 5: Tree graph for 



such expansions, logarithms of k III . Since we will use <3? in, „, „ 

^ ' ° ' n ' 9(k, P) of a parton target. 

for example, a calculation of hard scattering at scale Q, but 
need its value when \k 2 \ is much less than Q 2 , we must evaluate 



the right-hand side of Eq. ( 5.11 ) after a renormalization-group transformation using Eq. (5.2) 



Thus practical calculations will use 



+ o 



(5.12) 



On the right-hand side of this formula we have transformed the renormalization scale to 
\J\k 2 \. But the scale may be multiplied by a factor of order unity, which could be used to 
optimize for the likely size of uncalculated higher order corrections. 
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5.4.2 Integral of pCf 



The second result for the pCf was proved (with a slight variation) by one of us in [10|. It says 
that if we define a quantity / as the integral of <3? over k~ and kT, with an upper limit on kx 
and virtuality: 



f(x,M;g(/j,),m. 



J\k 



dk d kx 
(2tt) 5 



<f>(k,P;g(n),m,n), 



<\k 2 \,k^<M 2 

then this has an expansion in terms of the ordinary parton density: 



(5.13) 



f(x,M;g(/i),m,fj,)= / d££ T(£/x, M 2 //i 2 , g(//))/(£) + power-suppressed correction. 

J X 

(5.14) 

The formulation in [1C] differs only in that the limit was imposed only on kx and not on \k 2 \ 
as well, but the proof is unaltered. The lowest order value of the coefficient T is 5(£/x — 1), 



so that / = / + 0(g 2 ). From this follows the previously announced result Eq. (|5.6| ); the MS 
pdf is approximately the integral of an unintegrated parton density. 

This result gives a sum rule involving an integral over the whole non-perturbative region 
for the pCf, thereby providing a constraint on modeling of this region given that the ordinary 
pdf has already been measured. 

5.5 Recursive factorization for weighted pCf in terms of pCf 



In Eq. (4.37) the weighted DIS cross section was expressed with an integral over exact parton 
momenta in terms of the unweighted pCf. With the same reasoning, there is an analogous 
factorization property for the (weighted) pCf. Since we will wish to use it to give conditional 
probability densities for the internal partons in a MCEG, we normalize the formula to the 
unweighted pCf: 



*(k,P;ii) 



JV„=1 



dr + dr d A rr f dH$^ u {^, j) 



f(r+/P+;n) x 



^((r+,r~,r T ),P; M ) 
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f{f /p+]fi) 
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°°dM 2 J I (M 2 ;n>j,n)Q(M 2 /ni 2 J ) 
2tt Z{ji'j/n,m/n) 

dL(p Q ; P-r; m ph ) F(p Q , P; fj,) 



$(r,P;n) 



*((r+,r-,fr),P;M) 



$(r,P;/i) 
W(f). 



n 

3=1 



(5.15) 
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The main change compared with Eq. (|4,37| ) is that the differential hard scattering factor has 
been given a different name, dH§ instead of da. In addition, the symbol for the internal loop 
momentum has been changed from k to r, to avoid a conflict with the symbol for the external 
parton momentum. The differential hard scattering coefficient dH$ is defined in terms of 
subtracted hard scattering graphs in analogy to da in Eq. (f4.38|): 



dH*,N u 0, rfj) = dL(l u ; k + f,0) H<$,, Nu (p, j) Z(n'j/n, m/n) 



(5.16) 



The approximated internal parton momenta r and lj are defined as before, except for replacing 
q and k in the definitions in Sec. 4.3 by —k and r. 

Notice how we have used a different symbol //j for the factorization scale than we did 
in factorization for the cross section. When we employ Eq. ( 5.15| ), we will use the fact that 
the anomalous dimensions of the numerator and denominator on the left-hand side are the 
same, so that we can replace the value /i ~ Q used in factorization for the cross section by 

a value fi ~ \j\k~ 2 \ appropriate for the use of (5.15). We must also choose // j of the same 
order, so that the hard-scattering factor in this equation is perturbatively calculable, without 
large logarithms of k 2 / n' 2 j or of k 2 / fi 2 . 

It is also possible to write ( |5.15| ) in terms of unapproximated momenta, as we did for the 



cross section in Eq. ( 4.43 ) 



®[w](k,P;yt) 



oo N u 

En 

N u =lj=l 
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x W(f), 



d% 
(2^)6 



Jl{ll 



H<s>,n u (ji, $(r, P; n) A Nu (l u ,r) 



dL{P ] P-r; m ph ) F(p Q , P; /i) 
$(r,P;n) 



N u 

n 

3=1 



(5.17) 



with r = Y^j lj + k. 

The lowest-order value of H^^i is the essentially the same as the coefficient for the pCf 
in terms of integrated parton densities: 



Hq>i(k,f) 



2irg 2 
\k 2 \z 



5(-l + 2k-r/k 2 ) +0(g 4 



(5.18) 



The lowest-order value of the H$ n u with more external outgoing partons is of order g 2Nu . 



6. Implementation 



We now put the information together to give an algorithm for a MCEG. First, the ordinary 
factorization formula ( |4.l| ) gives us a calculation of the inclusive cross section (and hence of 
the structure function) as a function of x and Q. Then ( 4.37 ) expresses the weighted cross 
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section in terms of the pCf. Exactly similar factorization formulas apply to the pCf: Eqs. (5.8) 
and ( |5.15| ). Lowest order values for the coefficients are presented below the corresponding 
factorization formulas. 

We now assume that we have started to generate an event by generating values of x and 
Q according to probabilities given by the usual factorization formula for the inclusive cross 
section. We wish to generate a full event. 

1. Given the momenta of the photon and the hadron, q 11 and P^, we generate a value 
N u for the number of final-state partons and then a final-state partonic configuration 
(ll, . . . , In u ) with the density extracted from Eq. (f4.43| ) 

n 



(2vr) e 



H Nv *(k,P) 

■A Nu . (6.1) 



a[l] 



j-- 

This can either be done directly, or by separate steps of generating a massless config- 
uration (h, ■ ■ ■ , InJ) together with values of k~ , hr and Mi, ■ ■ ■ Mjy u , as at the end of 
Sec. ^, with probability densities extracted from from Eq. ( 4.37| ). In the second case we 
need to calculate the true partonic variables. 

The hard scattering is evaluated perturbatively with the renormalization and factoriza- 
tion scales fj, and fij chosen to be of order Q, to avoid large logarithms. 

2. For each final-state parton j that is above some threshold of perturbative invariant mass 
lj > Qq, generate a final-state parton shower using the algorithm given in Q in its own 
CM frame. Then boost it to the true lj. Since we already know the exact momentum 
of the parton, we should omit the step of finding the invariant mass. 

Otherwise the parton is left unchanged for nonperturbative hadronization. 

3. If the initial-state parton is above the threshold of perturbative invariant mass \k 2 \ > Qq, 
generate its parton shower. This is done by looping back to the first step, with the 
photon momentum q replaced by the parton momentum —k, and with the factorization 
formulas changed to the ones for the pCf. 

Otherwise the parton is left unchanged for nonperturbative hadronization. 

4. Once all remaining partons that are above the perturbative threshold have been show- 
ered, apply nonperturbative hadronization. 

This recursive algorithm for event generation reduces the whole problem to a repeated 
application of the same procedure, so that the program size is independent of the number of 
particles generated and the computational resources are linear in the number of particles. 

7. Zero- and one-loop corrections to coefficient functions 

In this section we give the coefficient functions for the structure function and the parton 
correlation function up to one-loop order. 
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7.1 Structure function 



For the hard scattering coefficients for the DIS structure function we could work directly with 
diagrams for the hard-scattering component of factorization, i.e., diagrams for the subgraph 
labeled "H" in a decomposition like Fig. |||. However, we will find it convenient to put the 
graphs in the context of their use in the factorization theorem. Thus we consider graphs of 
the form of Fig. ^, where we have attached a full amplitude for the pCf $ to the graphs from 
which the one-loop hard-scattering coefficient function is calculated. 
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Figure 6: One-loop graphs for DIS. 

There are certain complications in the interpretation of these graphs. Observe first that 
there is some apparent double counting. For example, the ladder graph (b) also implicitly 
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Figure 7: A graph like this with an initial- 
state parton self energy is not used. 



Q 













Figure 8: In a graph like this the final-state 
self energy contributes to the hadronization 
of the final-state parton, and such graphs are 
disregarded in computing the hard-scattering 
coefficient. 




Figure 9: Approximation for region Ri of 
Fig. §(a). The small square brackets indi- 
cate where the collinear approximation of Eq. 



(4.14) is made. 



appears in (a), with the rung I2 of (b) being inside the $ subgraph of (a). We actually mean 
that we consider each term in Fig. || or |6] not as a graph as such but as a decomposition of a 
graph for the purposes of examining its behavior: a graph can have multiple decompositions. 
We will indeed sum over the different decompositions, and the subtraction procedure in the 
hard scattering will take care of the double counting. 

Next, whereas we have used an explicit graphical factor <J? to represent the connection to 
the target proton, we have omitted corresponding fragmentation factors for the connection of 
the final-state partons to the final-state hadrons. A more thorough presentation would have 
done this, but our focus in this paper is on issues associated with initial-state partons. 

Also, observe that self-energy corrections on the incoming parton line are always included 
in the pCf. Therefore we omit from our considerations decompositions such as Fig. [?]; the 
self-energy on the incoming line does not contribute to our calculation of the one-loop hard- 
scattering coefficients. 

Similarly self-energy corrections on final-state parton lines, as in Fig. [8L are always to be 
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considered as part of the hadronization, and these self-energies make no direct contribution 
are therefore disregarded in the calculation of the hard-scattering coefficients. But observe 
that in the application to an algorithm for a MC event generator, there are normalization 
factors associated with final-state parton in the definition of the properly normalization hard 
scattering cross section d&N u - These are the factors of Z in Eq. ( f4.38| ), and their calculation, 
given in [Q], does indeed include a consideration of self-energy graphs. 



7.1.1 Zero- loop graph (a) 

Graph (a) has a tree graph coupled to the full pCf. The hard scattering contribution corre- 
sponds to the region where the incoming parton line k is collinear to the target. Comparing 
the graph with the general form of factorization, ( |4.31| ) and (4.37), we see that the hard 
scattering coefficient is 

m = e 2 [1 + 0(g% (7.1) 

with a corresponding cross section 

do x = —5(k+ - q+) Ke 2 [1 + 0(g 2 )} = ^-5(k+ - xP+) Ke 2 [1 + 0(g 2 )}. (7.2) 

7.1.2 Graph (b) 

Now we treat the ladder graph, (b) of Fig. |(| We always suppose k is collinear to the target. 
Then there are two leading regions. In region Ri, fci, I2 and k are almost collinear to the 
target, with m 2 ~ \k 2 \ < \kf\ <C Q 2 . In region R2, k\ and I2 are at wide angle, with 12 
m 2 ~ \k 2 \ < \k 2 \ < Q 2 . 

First we must remember that, as noted earlier, graph (b) also appears as a particular 
case of graph (a). We constructed an approximation to graph (a), and obtained the lowest- 
order hard scattering. In terms of graph (b), this approximation corresponds to the case 
that momentum k\ is almost collinear to the target, i.e., it is a good approximation in region 
R\. Therefore our analysis of the regions for (b) will be designed to obtain an appropriately 
subtracted hard scattering corresponding to the other region. 

This is done in terms of the approximation operations for the two regions of graph (b), 
which correspond to different projections of the parton momenta on massless momenta. The 
projection for region R\ we denote by 13 (k\, Zi, I2, k) with Jacobian Ax, and the projection 



12 Here, for convenience of power counting, we choose m to be the typical small scale. In general all we need 
is some scale s which is much smaller than the hard scattering scale Q. That means, the virtuality of the 
external partons of the hard subgraph is of order s 2 . If s is much larger than the hadronization scale A then 
we can recursively apply factorization to the jets generated by those external partons. 

13 Notice that no approximation is made on k or I2 in this region; these momenta are considered as part of 
the pCf. 
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for region R2 we denote by (ki, h, h, k) with Jacobian A2. The values are: 

k = (-g + ,0,0 r ) , 

h=(0,q-,d T y (7.3) 



and 



Ax = 1 + K/q- = 1 + 2xP + k^/Q 2 , (7.4) 
('it — hr) 2 Q~ + h ~ '2 'it ~~ hr \ 



r _ I ('it - hr) 2 <T + '2" ~ fT fgr ~ 'it \ 
2 " l v 4(g- + / 2 --/ 1 -)' 2 ' 2 )> 

~ I , (Jrr - hr) 2 —q~ + h - h 'it ~ hr \ 

h -{- q + 4 (<1 - + ir-i 2 -)' 2 '^2— J- 

Before approximation, graph (b) is represented by: 



(2vr) 6 ^ 6 )(/ 2 -A ; + A :i ). (7.7) 



d% 
(2tt) 5 2E 2 

Our goal is to apply the general decomposition by regions that we gave in Eqs. (|4.4j ), ( [4.25 ) 
with our chosen definition of the approximation operator Tr. We will construct terms T b \ = 
Cm(T b ) and T b 2 = Cjn(Tb) and will explicitly check that the sum is a good approximation 
to the original graph: 

T b = T bl + T b2 + 0(\T b \m 2 /Q 2 ) (7.8) 
everywhere that k is collinear to P. 

Region R\. In this region, \k\\ is 0(m 2 ), so I2 and k are almost parallel. The corresponding 
approximation, symbolized in Fig. ^, is 

dk+ r d 5 k 



T w = J T m T b = J— J __$(i)(^,P) dL(l 1 ;q + h)H l (i 1 ,ki;q,m = 0), (7.9) 
h 

^Hh,P) = J j^eHk,P) dL(h;k-k l) {kl _ 9 m2)2 (7.10) 
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representing a one-rung iteration of the pCf. The approximated momenta are defined by Eq. 
( |7.3|) . The result in Eq. ([O]) gives of course just an example of a lowest order hard scattering, 
and we already know that the contribution to the hard scattering factor is Hi(li, fei; q, m = 
0) = e 2 . 

When constructing the subtractions for region R2, we need the formula for T b i in terms 
of the original phase space integration: 

/r fj&h r f ]6i.i p 2 2 

rw = J wr m dm kl + q) dL{h] k ~ kl) Al Wi -m 2 ) 2 • (7 ' n) 

Also we need to verify that T b i is indeed a good approximation to T b in region R\\ 

/f d e k f d e ki 

(r„-r M ) = J — ^#(*,P) J j^ g dL(h;k 1 +q)dL(l 2 ;k-k 1 ) 

X 5T^ (1 - Al) ' <7 ' 12) 

Since 1 — Ai = k^ /q~ , which is small in R\, the desired suppression follows. 

Region R2. In this region, both of the final-state partons l\ and I2 are at wide angle. The 
approximation for the whole graph in this region is 

/■ ^ f dk+dk-d 4 k T „, , r .~ - ~ 
jT R2 T b = J ( 27r )6 ^(fc,P)d^(/i,/2;g + fc;m = 0) 

_^5W(Z a _(Ai-fc))F 2 . (7.13) 



Here #2 = e 2 g 2 / (kf) 2 , while fc, fei, Zj, and Z2 are given by (7.5). 

The hard scattering factor H2 has a collinear singularity when I2 and k become parallel. 
However, this is in the smaller region R±, where T b i is a good approximation, and for which 



we are required to apply a subtraction, according to the general formula ( 4.25 ). Therefore we 
define to be the approximation of the original graph T b minus the approximation in R±, 
i.e., Tfei. Then the singularity in the collinear region will be canceled point-by-point in the 
integration space: 

J T b2 = J V{-k 2 /m 2 )T R 2(T b -T bl ) 

- 1 wF* {k - p) I Br* I (OT (2 ^ (6,(!l - q - kl) (7 ' 14) 



*b i2K)9sW{h _ k + ki) v[ _~ k 2/ m 2 ) L* 

(27r) 5 2£ 2 (k 2 ) 
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Here A' lb is the Jacobian in the subtraction term with the momenta k, etc. replaced by the 
approximated momenta k\, etc., that are appropriate for region R 2 , i.e., 



A' 



ib 



Ai 



q - k~/2 

Zf - k-/2 + l 2 ~ 



q - k~/2 



2(if - k{ 



Ai 



'i 

, m 2 

1 + 



o 5 



(7.15) 



The Jacobian A2 is given by Eq. (4.22) with N u = 2. The collinear veto factor V(ki/m), 
defined in Eq. ( [4.29| ), is used to remove a spurious nonleading collinear divergence in the 
differential cross section due to the factor H 2 = e 2 g 2 /(k 2 ) 2 . Since T b2 is designed to handle 
momentum configurations far from the collinear region, the precise functional form of the 
veto factor in the collinear region is irrelevant to physics. Any other function which is zero 
when \kf \ <C m 2 and unity when \kf \ 3> m 2 would be equally good. 

It follows from Eq. (7.14) that the contribution of the one-rung ladder to the hard scat- 
tering is 



V(-k 2 /m 2 )^(l-A' lb ) 



(*?)' 



(7.16) 



Validity of approximations We can check that is indeed power suppressed in region 



r b2 = VT R2 (r b -r bl ) = o |r 6 -r, 



61 



(k 



2\2 



2\2 



where we have used Eq. ( |7. 12| ) to estimate the size of T b — T b i. In region R 2 , where \k 2 \ is 
the typical scale that we ignored and \k 2 \ is the typical hard scale, we can see the following 
is true. 

\k 2 



b2 



(T b -T bl ) = 



k\ 



bi 



(7.18) 



Combining the error estimates in different regions we can see that T^i and 1^2 indeed 
satisfy our requirement that = F b i + T b2 + O^T^m 2 /Q 2 ) in all regions. 



7.1.3 Graph (c) 

Next we treat graph (c) of Fig. ||. The treatment is essentially the same as that of graph (b), 
except that the subtraction is about final-state fragmentation. The graph has two leading 
regions. For its smallest region, for which we use the notation R&, p± and p 2 are almost 
parallel. In region R c2 , pi and p 2 are at wide angle. Since final state subtractions have been 
treated in detail in we will give the result directly without details. The major difference 
is in the projections to massless momenta in DIS and e + e~. Again, we use symbols like l\ for 
approximated momenta in region R c \ and symbols like l\ for those in region R c2 . 
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The graph can be represented as 



eV 

9 (7.19) 



(q + k) 2 — m 2 

Region R&. In region R c \ of graph (c), T c is approximated as 
j = J T Rcl T c 

= IT I HT^*^ f -^K(2^\h-q-k) 
J 2vr J (2vr) & J ( 27r )5 2 |Z 1 | 

fdM 2 f_fpi f d 5 p 2 6 (6) e 2 g 2 9(M 1 2 / / x 2 J ) 

x y (2^/ (^F^ (27r) ^ ()( ^- pi - p2) [( P1+ p 2 ) 2 -m 2 ] 2 - 

e 2 g 2 e(M 2 /fi 2 )A cl 



(q + /c) 2 — m 2 



(7.20) 



where the Jacobian is A c i = — A; - ) and the mass of parton ?i is M 2 = if = {p\ + p 2 ) 2 - 

This value is in fact a full pCf times the lowest-order hard scattering times an order g 2 
approximation to final-state showering. 14 Thus it is to be found as a contribution to full 
factorized cross section with the lowest-order hard scattering; there is no new contribution to 
the hard scattering. 

Region R c2 . For region R c2 of graph (c), we have the approximation 

eV 

9 (7.21) 



[(Pi+p 2 ) 2 ] 2 ' 

which is divergent in the collinear region R c \. After subtraction we have 
Jr c2 = J V(M 2 /m 2 )T Rc2 (T c -T cl ) 



(2tt) 6 

[(Pi +P2) 



> ] ^llm 2 ) u J 2g2 \\ - e((pi + p 2 ) 2 /^)A' lc ] A c2 , (7.22) 



14 Note that, whereas we drew graph (a) with an explicit factor for the initial-state interactions, we did not 
explicitly give it its corresponding final-state factor. 
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where A' lc is the Jacobian calculated with massless momenta for region i? 2 

A' lc = J}t P \ = P l^~ol = + 0{k-/lT)) = A 1 (1 + 0(m 2 /Q 2 )). (7.23) 
p 1 +p 2 - k~ p 1 + p 2 - 2k 

A c2 is defined using Eq. (|4.22Q with N u = 2. It is easy to check that T c = T cl +T c2 +0(m 2 /Q 2 ) 
in all regions. 

It follows that there is the following contribution to the hard scattering H 2 : 



^iV) r , e ,^ 212 



+V2Y 



(Pi + P2 



(7.24) 



7.1.4 Graphs (d) and (e) 



Graphs (d) and (e) of Fig. [6| are hermitian conjugates of each other. They each have only one 
leading region, where the two final-state particles are at wide angle. We have 

fr<= I As*M (Pvli # / 1 # (2*?^(ii - q -h) 

J J (2vr)6 V J J (2n) 6 J QvfZEx J (2^f2E 2 K ' yi q l> 

x^Wfo ~ k + h) {k2 _ m2 ^f +k)2 _ m2] (7.25) 

J T e = H.c. Jr d = J r d , (7.26) 
/" /" f d 6 k /7 , f d 5 h f d 5 l 2 



(2vr) 5 2£i J (2tt) 5 2E 2 
e 2 g 2 A 2 

(g + fc)2(Z* 2 -fc)5 



x (2vr) 6 ,5( 6 )(/ 1 + Z 2 - 9 - k)- f^f 2 (7.27) 



where A 2 and the approximated momenta (k, l 2 , • • • ) are defined using Eq. ( |jg) and Eq. 
(pUSD with iV M = 2. 

The contribution to the hard scattering of these two graphs is 

eV 

2 4 — . (7.28) 

(q + k) 2 (l 2 -k) 2 

7.1.5 Virtual corrections: graphs (f) and (g) 

Graph (f) and (g) in Fig. || provide vertex corrections to one-parton production. They have 
only one leading region, where the loop momentum is hard. There is also a UV divergence in 
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the loop integral which we renormalize using the MS-scheme, so that 
d 6 k „. f d% 



+ UV counterterm, 



(2tt) 6 (fcf - m 2 ) [(g + fci) 2 - m 2 ] [(A; - fei) 2 - m 2 ] 
r 9 = H.c. y T g (7.30) 



n 



d 6 fci eV 



+ MS counterterm, 



(2vr) 6 kKq + hyik-h) 



(2tt) 6 v ' ; 7 (2vr) 5 £i 

xAi-^fln^-sV (7.31) 
128vr 3 \ n 2 ) y ' 

Note that the vertex correction depends on the renormalization scale /u; this just shows that 
the "current" j = l/2[<^> 2 ] has an anomalous dimension. The renormalization group of the 
current can be found in . The anomalous dimension is canceled by the anomalous dimension 
of e 2 . 

7.1.6 Total 

Combining the 0(g 2 ) hard-scattering coefficients from graphs (b) to (e) in Fig. [f| we then 
have the hard-scattering coefficient for 7* + parton — ► p\ + p 2 - 

H 2 = V(- (fc-p 2 ) 2 /m 2 ) r/r eV oi - (1 - A' 16 ) + 2< ' V 



[(fc-p 2 ) 2 ] 2 i0/ (g + fc) 2 (p2-fc) 2 

+ V{{p l+ p 2 f/m 2 ) [1 - 6((pi +p 2 ) 2 A4)A' lc ] + O(eV), (7.32) 

[(Pi +P2YV 

with A' lb and A' lc defined in Eq. ( [7.15 ) and ( 7.23j ), respectively, and V being the veto function 
that cuts out the singularities. The approximated momenta pi, etc. are defined in Sec. 4.3 
with N u = 2. 

From graphs (a), (f) and (g) we find the one-parton term with one-loop corrections: 



Hi = e 2 



7.2 Parton correlation function 

In our model theory, the current for the "photon" is just a scalar vertex with two parton 
lines. Graphs for the parton correlation function are the same as for the structure function 
multiplied by: 



- 41 - 



2 

• Two full parton propagators at momentum k, i.e., a factor l/(k 2 ) at leading order. 

• A factor of g 2 instead of the factor e 2 associated with the vertex for e 2 [0 2 ]/2. 

• A factor —2ir/k 2 to remove the prefactor Q 2 /(2ir) in our definition of the structure 
function — see Eq. fl2.4|). 

8. Conclusions and future work 

We have shown that in a full treatment of a Monte-Carlo event generator the derivation 
of standard factorization is not sufficient. The physics in standard factorization is indeed 
unchanged. What is at issue is the validity of the identification, for example, of the momenta 
of the external partons of a lowest-order hard scattering with the naive parton model values, 
Eq. (2.5). It is not just that showering of the parton k in Fig. |3| gives it nonzero values k~ 



and kx, for that is implicit already in a correct definition of a parton density. The problem 
that forces us to change views is that k + is changed by the showering of the other parton. 
We have shown how these issues can be handled by the use of a new factorization theorem 
that uses unintegrated parton correlation functions, and in a model theory we have derived 
an appropriate theorem. 

In the new factorization theorem, we consider the structure function differential in the 
full final state, or more conveniently we consider the structure function with an arbitrary 
weight function applied to the final state. The form of the factorization enables a systematic 
first level of decomposition of the partonic structure to be made, differential in exact parton 
kinematics. The hard-scattering has a systematic perturbative expansion to all orders of 
perturbation theory. Instead of ordinary parton densities, the theorem uses parton correlation 
functions. 

Factorization also applies to the parton correlation functions, thus enabling a system- 
atic recursive decomposition of the fully exclusive parton structure of generated events. The 
formalism covers arbitrarily nonleading logarithms with the use of an expansion of the coef- 
ficient functions in powers of the coupling, without logarithms, and with the use of suitable 
renormalization-group transformations. It lends itself very naturally to a recursive implemen- 
tation suitable for a MCEG. 

The new factorization theorems that we have just summarized give a decomposition of 
the final-state into its exclusive components. We also need separate factorization theorems 
that give the cross section and the parton correlation function integrated over final states; 
these enable the first pair of factorization theorems to provide conditional probabilities to 
be used for the MCEG. For the cross section integrated over all final states, factorization is 
just ordinary factorization, which expresses it in terms of ordinary parton densities. We then 
showed that an almost identical argument gives a very similar result for the parton correlation 
function. 

The formalism is completed with the aid of renormalization group and DGLAP evolution, 
that allow perturbative calculations to be done at the appropriate scales in each component 
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of the different factorizations. The same ideas apply to final-state showering, but in a simpler 
way g]. 

Using the subtractive procedure proposed in || to define higher order corrections to 
the kernels, we constructed an algorithm for a MCEG for DIS (including initial-state parton 
showering), in which coefficient functions for both the hard scattering and the showering can 
be obtained at arbitrarily non-leading order. The subtractions are applied point-by-point in 
momentum space, so that the hard scattering coefficients and the evolution kernels are normal 
functions, unlike the singular generalized functions that appear in inclusive calculations from 
conventional factorization. 

Because of the subtractions, the hard scattering coefficients are not necessarily positive. 
An EG with weights is in general to be expected Jjj]. An algorithm for weighted EG was 
proposed in [12|. 

Although our derivations only apply to non-gauge theories, and there is much work to be 
done to extend our work to full QCD, we believe our work indicates how one should approach 
the full QCD problem. One considers the regions of momentum space that dominate but 
without relying on any unitarity cancellation of a soft region. This is like the situation in 
factorization with transverse-momentum-dependent (TMD) parton densities |11, 12], which 
work applies to inclusive cross sections that are sensitive to parton transverse momentum. 
In QCD (and any gauge theory), two complications arise. The first is that there is an extra 
non-canceling soft factor to be treated. The second is that to be gauge- invariant, the parton 
correlation functions must be defined with Wilson lines in the operators, and factorization 
requires [13, 14, 15, 16, 17] the Wilson lines to be in certain directions. Examination of 
the general techniques used in these works suggests that they should continue to apply to 
the weighted cross sections we use to understand the final state. Naturally the appropriate 
factorization must involve gauge-invariant unintegrated parton correlation functions, so that 
parton momenta can be treated exactly. There are extra techniques already available for 
handling the extra complications, for example the Collins-Soper evolution equation for the 
direction dependence of the Wilson lines. So many of the elements needed to extend our 
results on MCEGs to full QCD are available, and simply need to be adjusted and applied in 
the new context. 



Recently Collins and Metz [17] have performed a detailed general analysis of the directions 
of the Wilson lines needed for factorization; technically it depends on an understanding of the 
allowed directions of certain contour deformations in virtual graphs that are needed in the 
proof of factorization. In the case of both e + e~-annihilation and of DIS, their analysis of TMD 
factorization does not have any requirement that the cross section be inclusive. Since such 
cross sections involve a soft factor, no unitary cancellations are needed to get factorization. 
But one will probably need multi-parton soft factors to accompany multi-parton production 
processes, with corresponding new nonperturbative features. 

It is the use of TMD factorization that saves the need for the unitary cancellation of the 
soft factor that happens in non-TMD factorization. 

But as Collins and Metz explain, TMD factorization in hadron-hadron processes, like 
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Drell-Yan, continues very critically to need a unitarity cancellation. Therefore, the limits 
to MCEGs are likely to be more severe here. This suggests that a full treatment will need 
a discussion of spectator-spectator interactions that goes far beyond a normal factorization 
framework. Elements of this are undoubtedly already present in some existing EGs like 
PYTHIA. 
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